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Abstract  -  The  performance  of  a  tracking  filter  can  be 
evaluated  in  terms  of  the  filter’s  optimality  conditions. 
Testing  for  optimality  is  necessary  because  the  estimation 
error  covariance  as  provided  by  the  filter  is  not  a  reliable 
indicator  of  performance,  which  is  known  to  be 
“optimistic”  ( inconsistent )  particularly  when  there  are 
model  mismatches  and  target  maneuvers.  The 
conventional  root-mean  square  (RMS)  error  value  and  its 
variants  are  widely  used  for  performance  evaluation  in 
simulation  and  testing  but  it  is  not  feasible  for  real-time 
operations  where  the  ground  truth  is  hardly  available. 
One  approach  for  real-time  reliability  assessment  is 
optimality  self  online  monitoring  (OSOM)  investigated  in 
this  paper.  Statistical  tests  for  optimality  conditions  are 
formulated.  Simulation  examples  are  presented  to 
illustrate  their  possible  use  in  evaluation  and  adaptation. 
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1  Introduction 

From  design,  procurement,  to  deployment,  performance 
evaluation  remains  one  of  the  issues  that  target  tracking  and 
automatic  target  recognition  (ATR)  system  engineers, 
government  buyers,  and  systems  operators  must  address  in 
different  phases  of  the  systems’  life  cycle.  In  design  and  test 
stages,  performance  evaluation  is  relatively  easier  not  only 
because  the  ground  truth  is  available,  but  also  because  the  tests 
can  be  repeated.  During  deployment,  however,  the  ground  truth 
about  a  target  can  never  be  known  a  priori  and  there  seldom 
exists  a  second  chance  in  a  time-critical  hostile  environment  to 
verify  algorithm  performance. 

There  are  at  least  two  types  of  performance  evaluation.  One  is 
comparative  (and  competitive)  among  a  variety  of  algorithms  or 
designs  wherein  a  set  of  performance  metrics  are  used  to  rank  the 
algorithms  or  designs  for  one-out-of-many  selection.  It  typically 
requires  the  ground  truth.  In  most  cases,  no  single  algorithm  or 
design  would  prevail  for  all  performance  metrics  and  no  single 
metric  can  capture  all  performance  aspects  under  all  conditions 
(i.e.,  scenarios).  The  other  type  of  performance  evaluation  is 
absolute,  which  is  concerned  with  a  particular  algorithm  or 
design  in  run  time.  The  second  type  of  performance  evaluation  is 
the  focus  of  this  paper. 

Although  comparative  metrics  can  be  used  when  the  ground  truth 
is  available  [4,  18],  an  alternative  approach  is  to  make  use  of 
design  specs  and  optimality  conditions  as  a  means  for 
performance  evaluation.  Indeed,  modem  target  tracking  and  ATR 


system  designs  are  model-based  and  they  are  optimized  or  tuned 
with  respect  to  certain  design  criteria  or  cost  functions.  When  the 
design  conditions  are  met,  the  filter  is  expected  to  perform  within 
specs.  In  a  sense,  each  filter  ought  to  have  an  optimality  self 
online  monitoring  (OSOM)  module  to  ascertain  if  the  operation  is 
optimal  and  if  not.  to  provide  timely  feed  back  for  online  tuning 
and  adaptation.  In  a  sense,  this  forms  a  closed  loop  involving 
sensor  management,  thus  related  to  such  research  areas  as 
adaptive  control  and  fault  detection  and  isolation  (FDI). 

This  paper  presents  an  OSOM  system  for  target  tracking  filters 
and  their  fusion  via  adaptive  sensor  management.  It  monitors  the 
optimality  conditions  of  a  tracking  filter  using  statistical  tests 
(e.g.,  normality  and  whiteness)  based  on  the  innovation  sequence 
and  relates  the  deviation  from  the  optimal  conditions  to  the  errors 
developed  in  the  target  state.  It  can  further  calculate  the  actual 
process  and  measurement  noise  covariance  matrices  as  compared 
to  those  used  by  the  filter.  Together  with  the  sensitivity  of 
optimality  conditions  with  respect  to  various  filter  parameters 
(e.g.,  Q,  R,  P0 ),  it  enables  on-line  tuning  or  adaptation  so  as  to 
maintain  the  optimality  conditions. 

The  paper  is  organized  as  follows.  In  Section  2,  a  typical  tracking 
filter  is  described  together  with  analytic  performance  prediction. 
In  Section  3,  a  series  of  optimality  tests  are  presented  for 
performance  monitoring.  In  conjunction  with  OSOM,  methods 
for  filter  tuning,  adaptive  distributed  fusion,  and  active  sensor 
management  are  outlined  in  Section  4.  Section  5  concludes  the 
paper  with  some  future  work.  Simulation  examples  are  presented 
to  illustrate  the  OSOM  procedure  in  the  paper. 

2  Tracking  Filter  &  Performance 

Tracking  filter  design  strikes  a  balance  between  noise 
performance  and  dynamic  responsiveness  to  target  maneuvers. 
Given  the  sensor  error  characteristics,  the  design  boils  down  to 
the  choice  of  covariance  for  process  noise,  which  is  a  modeling 
tool  widely  used  in  the  Kalman  filter  to  deal  with  uncertainty. 
However,  the  estimation  error  covariance  of  the  Kalman  filter 
becomes  inconsistent  when  large  maneuvers  occur.  The  sensor 
noise  only  covariance  and  steady-state  biases  or  lags  in  position 
and  velocity  may  be  used  instead.  Design  optimization  then 
consists  of  finding  the  filter  gain  given  the  maximum  possible 
maneuver  that  minimize  the  overall  root  mean  square  (RMS) 
errors  [6]. 

2.1  Target  &  Measurement  Truth  Models 

Consider  the  following  target  model: 

x4+i=Fixi+Gl(uI+'T)  (la> 

z*=H*xt+wt  (lb) 
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where  for  the  one-dimensional  (ID)  case  (i.e.,  in  one  coordinate, 
say,  the  x-axis)  the  following  variables  and  matrices  are  defined: 


H4  =  [l  0] 


(2a) 

(2b) 

(2c) 

(2d) 


T  is  the  sampling  interval,  Ua  is  an  unknown  deterministic 
acceleration  input,  \k  is  a  random  acceleration  process  with  zero 
mean  and  variance  crv2,  zk  is  the  measurement,  and  wk  is  the 
measurement  noise  with  zero  mean  and  variance  Rk  =  a„2. 


2.2  Kalman  Filter  &  a- P  Filter 


Given  the  target  motion  and  measurement  models  in  (1)  and  (2), 
the  corresponding  Kalman  filter  consists  of  the  time  update 
(propagation)  step  (if  uA.  is  known): 


=  FP„Fr  +Qt  =  FP<.,Fr  +  GkGka2 


(3a) 

(3b) 


where 


=  E{xk  I Z 


is  the  state  estimate  at  time  k  with 


measurements  up  to  k  denoted  by  Z*  =  {z0,  Z\,  ....  zk}. 


=F{(x4  -xtlt)(xt  -ittk)T  \Z 


is  the  estimation  error 


covariance. 


=  £{x,lIZi)  's  the  one-step  ahead  state 


prediction,  and  p  -  fii\  -i  itv  -i  \T  \7k\  is  the 

rk+Uk  ~  \*k+l  Aifc+1I*AA*+1  A-k+Uk)  '  ^  i 

prediction  error  covariance,  and  the  measurement  update  step: 


*«*  =*«*-, +  Ki(zl-Hlitlt_1) 
PM=(I-K,Ht)P4l(_1 


(3c) 

(3d) 


where  the  Kalman  filter  gain  is  given  by: 


2.3  Analytic  Performance  &  Design  Analysis 

A  Kalman  filter  in  the  steady-state  (i.e.,  the  a-P  filter)  can  be 
viewed  as  a  linear  system  and  so  is  the  truth  model.  Together, 
they  are  driven  by  random  measurement  noise,  random 
acceleration  input,  and  deterministic  acceleration  (maneuver) 
with  unknown  in  time  and  magnitude.  When  Ua  =  0  and  \k  =  0, 
the  state  estimation  error  covariance  position  and  velocity 
components  contributed  to  by  the  sensor  noise  only  ( SNO )  are 
given  by: 


pp 

2p  -3aP  +  2a2  2 

(5a) 
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1 

(N 

1 

¥ 

p™ 
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(5b) 

a(A-2a-  P)T  " 

pr 

<N 

II 

(5c) 

a{A-2a-  P)T 

which  show  the  reduction  of  variances  of  the  filtered  state  as 
compared  to  the  variances  of  noisy  measurements.  Similarly,  the 
state  estimation  error  covariance  position  and  velocity 
components  contributed  to  by  a  constant  acceleration  of  level  A, 
when  uA  =  A  and  Va  =  0,  are  given  by: 

P*  = 

Pf  = 


T 2 

(1  -a)— A 
(^-0.5  )TA 


(6a) 

(6b) 


which  are  also  called  the  acceleration-induced  bias  or  lag  in 
position  and  velocity,  respectively. 

The  random  process  noise  is  typically  used  to  account  for 
modeling  error  and  uncertainty  about  target  maneuver.  When  the 
system  is  subject  to  a  random  acceleration  with  zero-mean  and 
variance  a„2,  the  state  estimation  error  covariance  position  and 
velocity  components  contributed  to  by  this  random  acceleration 
noise  {RAN),  when  Ua  =  0,  are  given  by: 


K,.=  Pai^H'S)1 

with  the  measurement  prediction  error  covariance: 
Sa  =H  P„  H)  +Rt 


(3e) 

(3f) 


In  the  steady  state,  the  above  Kalman  filter  for  the  motion  and 
measurement  models  (1)  and  (2)  becomes  the  so-called  a-/?  filter 
with  the  constant  gain  denoted  by: 


K 


k 


(4a) 


where  the  gains  are  given  by  [2,  3,  10,  16]: 

a  =  -0.125(772  +  8// -(77  +  4)-J  if  +8//)  (4b) 

P  =  2(2  -  a)  -  Wl  -a  (4c) 


pra, 

1  =  aal 

(7a) 

P™ 

'  =  —  a2r 

(7b) 

T 

P"“ 

,  P{a-0.5P)  2 

t2  n  w 

(7c) 

T2(\-a)  w 


which  are  the  solution  of  the  algebraic  Riccati  equation  (3b )  and 
(3d). 

The  gains  a  and  P  can  be  related  to  the  damping  ratio  £  and 
bandwidth  (undamped  natural  frequency)  co,,  for  a  continuous¬ 
time  filter  [10,  24]  by: 

a  =  2Cjp-p/2  (8a) 

P  =  (coT)2  (8b) 

0<or<l,  0<  P  <2(l-a)  (8c) 


using  the  tracking  index  (target  maneuverability  index): 


rj  =  T 2^~ 
cr 


(4d) 


which  can  also  be  viewed  as  signal  to  noise  ratio  (SNR),  i.e.,  the 
position  error  caused  by  constant  acceleration  T2<jv  vs.  the 
position  measurement  error  cjk.  In  the  steady  state,  the  transfer 

function  from  measurement  zk  to  state  estimate  *  is  also  the 

4  xaia 

Wiener  optimal  filter. 


The  last  condition  (8c)  ensures  the  filter  stability. 

Combining  (5)  and  (6)  gives  the  maximum  RMS  errors  in 
position  and  velocity  errors: 


RMS™'  =  +  Pp™ 

RMS™'  =  yjPT  +  Pfm" 

where  Amax  is  the  maximum  possible  acceleration  level. 


(9a) 

(9b) 
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2.4  Illustrative  Examples 

A  target  is  stationary  at  x,  =  0  from  t  =  0  to  200  sec.  It  then 
accelerates  at  u,  =  5  m/s2  from  t  =  200  to  400  sec  and  next  moves 
at  the  achieved  velocity  of  1000  m/s  from  t  =  400  to  600  sec.  It 
then  decelerates  at  u,  =  -5  m/s2  from  t  =  600  to  800  sec  and  finally 
stops  at  x,  =  400  km  up  to  t  =  1000  sec.  The  sampling  interval  is  T 
=  1  sec.  In  addition  to  the  deterministic  input  u„  a  random 
acceleration  v,  ~  iA/( 0,  0.12)  (i.e.,  <7r  =  0.1  m/s2)  is  also  used  in 
some  of  the  simulation  examples  presented  below. 

The  target  is  observed  in  position  by  a  sensor  with  its 
measurement  noise  standard  deviation  being  =  5  m.  The  initial 
estimation  error  covariance  is  P0  =  diag([  1002,  102]).  The  filter  is 
initialized  with  a  sample  drawn  from  JA/|0,  P0}. 

Fig.  1  shows  the  Kalman  filter  gain  as  a  function  of  time  vs.  a 
and  P  calculated  by  (4).  The  position  gain  is  slightly  slower  than 
the  velocity  gain  but  both  converge  to  the  steady  state  a  and  P 
within  50  seconds.  For  the  particular  design,  the  damping  ratio  is 
£  =  0.6398  and  the  natural  frequency  is  ft},  =  0.1345  rad/sec  or 
0.0214  Hz,  which  corresponds  roughly  to  a  time  constant  of  47 
sec  for  transient. 

Figs.  2  and  4  compare  the  RMS  errors  in  position  and  velocity 
calculated  from  50  Monte  Carlo  runs  (blue),  the  filter  provided 
error  covariance  (7)  (green),  and  the  predicted  maximum  errors 
(9)  (red),  respectively.  The  simulated  RMS  values  (blue)  match 
well  with  the  filter  predictions  (green)  in  quiescent  modes 
without  maneuver,  which  are  zoomed  in  Figs.  3  and  5.  During 
maneuver,  however,  the  filter  predictions  (green)  remain  flat  (i.e., 
erroneous  or  inconsistent  with  the  reality).  In  contrast,  the 
maximum  errors  predicted  from  (9)  (red)  with  =  u,  =  +5  m/s2 
match  well  with  the  RMS  errors  (blue)  over  the  maneuver 
periods. 

Figs.  6  and  7  show  position  and  velocity  errors,  respectively, 
where  a  single  run  (red)  is  superimposed  with  the  average  (blue). 
Typical  transient  behavior  for  a  2nd  order  system  is  visible.  The 
rising  time  is  about  45  sec,  consistent  with  the  system  bandwidth. 

Fig.  8  shows  the  innovation  sequence,  a  single  run  (red)  and  the 
average  (blue),  with  the  predicted  measurement  error  2cxbounds. 
The  blow-up  is  shown  in  Fig.  9  where  the  upper  and  lower 
bounds  are  in  dashed  green  and  red-colored,  respectively. 

When  v,  =  0,  the  target  is  driven  by  the  deterministic  ut.  The 
sensor  noise  only  error  predictions  (5)  are  shown  in  Figs.  10  and 
1 1  for  position  and  velocity,  respectively.  The  sensor  noise  only 
error  predictions  (5)  (dotted  green)  match  well  with  the  simulated 
RMS  values  (aqua)  whereas  the  filter  provided  error  covariance 
(7)  (green)  does  not,  which  is  pessimistic  (higher  than  the  actual). 
Fig.  12  shows  the  innovation  sequence,  which  is  well  within  its 
predicted  2  abounds  in  quiescent  modes  but  in  maneuver  periods. 

In  the  above  simulation,  the  analytic  formulas  match  well  the 
numerical  results,  indicating  their  usefulness  for  performance 
prediction.  However,  when  the  actual  value  of  crv,  <JW,  and  Amax 
differs  from  the  one  assumed  in  the  filter,  the  actual  tracking 
errors  will  be  different  from  those  given  by  the  formulas. 

This  is  shown  in  Figs.  13  through  15  for  the  position  error, 
velocity  error,  and  innovation,  respectively,  where  the  filter  uses 
( o;,) modd  =  2x(o;,)actua|.  As  shown,  the  predicted  values  and 
bounds  are  larger  than  the  actual  ones.  The  filter  is  pessimistic  or 
conservative.  The  opposite  is  also  true  when  the  filter  uses  values 
smaller  than  the  actual  (optimistic).  In  other  words,  the  prediction 
will  fail  when  there  is  a  mismatch  in  model  parameters.  Hence,  it 


is  necessary  for  a  filter  to  verify  its  design  parameters  in  run  time 
so  as  to  ensure  the  operational  optimality. 

3  Online  Monitoring  of  Optimality 

The  performance  of  a  tracking  filter  is  expressed  in  terms  of 
various  metrics.  There  metrics  are  evaluated  over  an  observation 
interval,  characterizing  how  far  away  the  state  estimate  is  from 
the  true  state  in  one  way  or  another.  Most  of  these  metrics  need  to 
know  the  true  state  for  their  calculation.  Such  metrics  are  good 
for  design  simulation  and  test  experiments  where  the  true  state  is 
available. 

In  this  paper,  we  are  interested  in  those  metrics  that  can  be 
implemented  by  a  tracking  filter  without  knowing  the  ground 
truth  [19].  One  approach  is  to  explore  the  filter  optimality.  A 
well-designed  (tuned)  Kalman  filter  should  provide  an  optimal 
estimate  of  the  state  (i.e.,  minimal  error  variance).  A  necessary 
and  sufficient  condition  for  a  Kalman  filter  to  be  optimal  is  that 
its  innovation  sequence  is  zero-mean  and  white. 

As  a  result,  these  properties  should  be  checked  routinely  to 
ensure  that  the  filter  is  operating  properly  [7].  If  so,  the 
estimation  error  covariance  can  be  trusted.  If  not,  further  tests  are 
needed  to  determine  what  may  go  wrong  (diagnosis)  and  on-line 
re-tuning  may  be  necessary. 

3.1  Tests  for  Zero-Mean  of  Innovation  Sequence 

We  first  consider  a  parametric  test.  When  the  innovation 
sequence  is  ergodic  and  Gaussian,  the  sample  mean  can  be  used 
to  estimate  the  population  mean  as  test  statistic.  Consider  an 
individual  component  of  the  innovation  denoted  by  e(k)  =  z{k)  - 

Z(k\k- 1). 

Denote  the  true  mean  and  covariance  of  the  innovation  by  m  and 
R.  And  denote  the  sample  mean  by  m  ,  which  is  estimated  by: 

1  k+N-l  D  /i  /~v\ 

m(k)  =  —  V  e(t)  ~  N{m,—} 

N  S  N 

where  N  is  the  number  of  samples.  Two  hypotheses  can  be 
formulated  as: 


Ha  m  =  0  (11) 

//,  m  i1  0 


The  probability  of  rejecting  the  null  hypothesis  Hn  at  the  a 
significant  level  with  threshold  ris  given  by: 


n 


Pr/n 


T—m  , 

>  ,  1  =  a 


Pr/n 


(12) 


The  zero-mean  test  on  each  component  of  the  innovation  is 
formulated  as: 


„  \>T  Reject  H0  (13) 

m  \ 

[<  r  Accept  //, 

At  5%  significant  level  {a  -  0.05),  we  have: 

r  =  1.96.1—  (14) 

V  N 

where  the  sample  variance  is  estimated  under  the  null  hypothesis 
as: 


R(k)  = 


1 

N 


k+N-l 


2>2« 


t=k 


(15) 
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Fig.  1  Kalman  vs.  a-p  Filter  Gains 


Fig.  2  Monte  Carlo  Position  Errors 

qm,  rm,  qa>  rg  =  0.1,  5,  0.1,  5 


Fig.  3  Position  Errors  (Blowup  of  Fig.  2) 

qm,  rm,  qa,  r  =  0.1,  5,  0.1,  5 


Fig.  4  Monte  Carlo  Velocity  Errors 

qm,  rm,  qa,  r  =  0.1,  5,  0.1,  5 


Fig.  5  Velocity  Errors  (Blowup  of  Fig.  4) 

qm,  rm,  qg,  r  =  0.1,  5,  0.1,  5 


Fig.  6  Sample  Behavior  of  Position  Errors 

qm,  rm,  q  ,  r  =  0.1,  5, 0.1,  5 


Fig.  7  Sample  Behavior  of  Velocity  Errors 


qm.  qa,  ra=  0.1,  5,  0,  5 


Fig.  8  Sample  Behavior  of  Innovation 

qm.  qa.  ra  =  0.!,  5,  0,  5 


Fig.  9  Blowup  of  Fig.  8 

qm.  qa,  r  =  0.1,  5,  0, 5 


Fig.  10  Position  Errors  with  Sensor  Noise  Only  Fig.  1 1  Velocity  Errors  with  Sensor  Noise  Only 


Fig.  12  Innovation  with  Sensor  Noise  Only 


For  nonparametric  test ,  an  asymptotic  relative  efficiency  (ARE) 
[11]  efficacious  small-sample  (=  10)  robust  (insensitive  to 
outliers)  nonparametric  distribution-free  linear  rank  test  from  the 
Chemoff-Savage  class  [8]  for  testing  shift  in  location  is  the 
Mann-Whitney-Wilcoxon  (MWW)  two-sample  statistics  [11,  14]. 
Two  independent  random  batch  samples  n  and  m  are  taken 
consecutively  of  the  innovations  sequence,  representing 
populations  X  and  Y ,  respectively,  where  the  total  number  of 
samples  N  =  m  +  n,  the  distribution  model  is 
Fy(x)  =  Fx(x  -  0),  and  the  null  hypothesis  of  identical 


distributions  H0  :  0  =  0.  A  symmetric  version  of  the  MWW  [14] 
test  is  given  by: 


W  =  —  XXsgn,<  v,) 


(16) 


For  selected  small  sample  sizes,  the  cumulative  density  function 
( cdf)  of  W  is  tabulated.  Thus  a  given  value  can  be  stored  and 
compared  with  a  threshold  at  given  level  of  significance  to  accept 
or  reject  the  null  hypothesis.  While  the  test  requires  ranking  of 
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the  batched  samples,  for  small  sample  sizes  ( m  +  n  <  20)  the 
computation  time  and  delay  incurred  are  very  small. 

3.2  Measurement  Screening  and  Dynamics  Detection 


statistics.  The  user  may  reject  the  null  hypothesis  if  W  is  too 
small. 

3.4  Tests  for  Zero-Mean  White  Gaussian  Condition 


In  order  to  both  flag  and  eliminate  measurement  outliers  Tukey’s 
method  of  Hinges  [22]  can  be  used  to  interpolate  quartiles. 
Given  a  small,  moving  batch  of  ranked  data  samples,  the  median 
of  the  sample  is  computed.  The  data  are  divided  into  high  and 
low  groups,  and  the  middle  (median)  values  of  the  groups  are 
denoted  as  Q1  and  Q3,  the  first  and  third  quartile,  respectively. 
The  inter-quartile  range  (IQR)  is  computed,  i.e.,  IQR  =  Q3  -  Ql. 
Next  two  “fences”  are  computed  as  Fence(lower)  =  Ql  -  1.5xIQR 
and  Fence(upper)  =  Q3  +  1.5xIQR.  Values  outside  the  fences  are 
considered  outliers.  To  note,  this  test  was  used  with  great  success 
in  eliminating  sensor  caused  outliers  with  real-time  data  [15]. 

The  same  method  can  be  applied  to  the  innovations  sequence.  In 
this  case,  values  outside  the  fences  can  be  used  to  detect  dynamic 
changes  not  accounted  for  in  the  filter  model. 


A  simple  statistic  that  contains  all  the  innovation  information 
over  some  finite  window  of  length  N  is  the  weighted  sum  squared 
residual  (WSSR)  [7]  defined  as: 

p(k)=  Y.e(t)S-'(t)e(t)~Z2{Np}~W{Np,2Np}  (21) 

t=k_N+l 

where  p  is  the  dimension  of  e(fc)  and  S (k)  is  obtained  from  the 
Kalman  filter. 

Under  the  Gaussian  assumption  for  Np  >  30,  the  threshold  for  a 
level  of  significance  of  a  -  0.05  to  accept  the  null  hypothesis  is: 

r  =  Np  + 1 .96^1 2Np  (22) 

3.5  Skewness  &  Kurtosis 


3.3  Tests  for  Whiteness  of  Innovation  Sequence 

The  parametric  whiteness  test  is  performed  to  check  statistically 
that  the  innovation  sequence  is  a  white  sequence.  The  sample 
covariance  as  a  function  of  delays  or  lags  is  used  as  the  test 
statistic.  For  a  component  of  the  innovation  sequence,  the  sample 
covariance  function  is  given  by: 

1  k+N~[  ( 1  '7\ 

R(k’A)  =  ~N  '"(*(] [e(r  +  A)  —  m(k)\ 

where  A  =  0,  1,  ...,  is  the  delay  or  lag  of  the  covariance  function. 
The  normalized  covariance  test  statistic  is  defined  as: 


When  a  Kalman  filter  operates  correctly,  the  innovation  sequence 
is  white  Gaussian  of  zero  mean.  The  skewness  and  kurtosis  can 
be  used  to  characterize  the  deviation  of  the  data  sequence  from 
the  normality. 

Consider  a  random  variable  z.  Its  skewness  and  kurtosis  are 
defined  respectively  as: 

£{(;->».)-}  (23) 

3  E{(z-mz)2}112 

£{(z-in.)4}  (24) 

7i  E{(z-mz)2}2 


P(t  A) 


R{k,  A) 
R(k,  0) 


R(k,  A) 
R(k) 


7V{0,— } 
N 


(18) 


where  mz  =  E[z]  is  the  mean  value  of  z.  The  denominator  in  (23) 
and  (24)  is  the  variance  of  z  raised  in  power. 


For  N  >  30,  the  95%  confidence  interval  is  +i.96/V)v  •  For  zl  =  0, 
p(k,  Zl)  =  1  so  the  interval  is  1  +  1.96/ViV  •  For  Zl  ^  0,  p{k,  Zl)  =  0 
so  the  interval  is  +1.96 /V)v  •  Similar  tests  can  be  constructed  for 
the  cross-covariance  properties  between  innovation  components 
as  well. 

A  more  efficient  implementation  of  (16)  is  to  use  the  fast  Fourier 
transform  (FFT)  to  compute  the  covariance.  Furthermore,  the 
autocorrelation  function,  R(  t),  of  the  batch  sequence  can  be  used 
to  compute  the  related  power  spectral  density,  S(f),  in  order  to 
gain  an  insight  into  the  uniformity  (or  lack)  of  the  shapes  of  S(f). 

A  highly  sensitive  small  sample  nonparametric  test  for  whiteness 
is  the  Shapiro-Wilk  [8,  26]  test  with  the  null  hypothesis  that  a 
sample  X\,  ....  xn  came  from  a  normally  distributed  population. 
The  test  statistic  is: 

(Z"=ifl.-*w>2  (19) 

X‘=i  (A-*)3 


The  skewness  characterizes  the  degree  of  asymmetry  of  a 
distribution  around  its  mean  value.  A  positive  value  of  skewness 
corresponds  to  a  distribution  with  an  asymmetric  tail  extending 
on  the  right  of  the  mean  whereas  a  negative  value  indicates  an 
asymmetric  tail  to  the  left.  The  kurtosis  measures  the  relative 
peakednes  or  flatness  of  a  distribution. 

For  a  Gaussian  variable,  both  skewness  and  kurtosis  are 
identically  zero.  So  they  offer  a  measure  of  the  deviation  from 
Gaussianity. 

3.6  Sequential  Estimation  of  Noise  Characteristics 

First  and  second-order  moments  of  the  noise  processes  can  be 
estimated  based  on  the  state  estimates  produced  by  the  Kalman 
filter,  which  provides  a  means  to  check  on  the  filter  modeling  and 
its  proper  operation. 

Construct  Nx  samples  of  most  recent  state  estimate  errors: 
f  =  Y  — X  =  X  — K  X  (25) 

lk  k\k  a/:IA  I  Ak\k  k  A k  111  [ 


where  x(i)  (with  parentheses  enclosing  the  subscript  index  i)  is  the 
i'th  order  statistic,  i.e.,  the  ith-smallest  number  in  the  sample,  x  is 
the  sample  mean,  the  constants  are  given  by: 


[ti,  u„]  = 


nTV  1 

(rnYVniT 


(20) 


where  m  =  [mi,  ...,  m„  |T  with  its  elements  being  the  expected 
values  of  the  order  statistics  of  independent  and  identically- 
distributed  random  variables  sampled  from  the  standard  normal 
distribution,  and  V  is  the  covariance  matrix  of  those  order 


If  the  system  has  an  unknown  constant  forcing  function,  it  can  be 
estimated  as 


(26) 


Similarly,  a  bias  in  the  measurement  can  be  estimated  based  on 
Nz  most  recent  measurement  prediction  errors: 


N. 


Viz, 


+  G  4b  /)] 


(27) 
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Fig.  13  Position  Errors  with  Parameter  Mismatch  Fig.  14  Velocity  Errors  with  Parameter  Mismatch 


Fig.  15  Innovation  with  Parameter  Mismatch 


95%  confidence  test  for  zero  mean  of  innovation,  sliding  window  =  40 


test  for  normality,  sliding  window  =  40 


Fig.  16  Zero  Mean  Test 


Fig.  17  Skewness  &  Kurtosis 


quiescent  mode 


innovation,  m 


Fig.  1 8  Histogram  of  Innovation 


maneuvering  mode  95%  confidence  test  for  zero  mean,  white,  gaussian,  sliding  window  =  40 
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Fig.  19  Histogram  of  Innovation  in  Maneuver  Fig.  20  Zero  Mean  White  Gaussian  Test  x2 


95%  confidence  test  for  whiteness,  sliding  window  =  40  at  time  =  150 


Fig.  21  Normalized  Correlation  Function 


95%  confidence  test  for  zero  mean  of  innovation,  sliding  window  =  40 


test  for  normality,  sliding  window  =  40 


95%  confidence  test  for  whiteness,  sliding  window  =  40  at  time  =  150 


Fig.  22  Zero  Mean  Test  (Sensor  Noise  Only)  Fig.  23  Skewness  &  Kurtosis,  Sensor  Noise  Only 


Fig.  24  Correlation  (Sensor  Noise  Only) 


It  was  shown  that  the  bias  estimators  (26)  and  (27)  are  unbiased 
for  optimal  state  estimates  [23].  It  was  also  shown  that  the 
following  process  noise  covariance  estimator  is  unbiased  [12]: 

Qa  =  £ A-m  -  G A )(t/H  - G  A )T  ] 

-  1  VlF  P  ¥r  -P  1  (28) 

N  L  {tk-j^k-jbk-j  rk_j+ 1) 

ly  X  j= 1 

where 


(29a) 


1  N' 

Gfc-iUfc-i  =— —  X^-7+l 

-r  7=1 

The  measurement  noise  covariance  is  estimated  by: 

1  A 

Ra  =tt— tZ[(Za-a«  AX  Am  -rkY 

iyz  1  7=1 

-  ^ (Ft_A_,F Ij  +  Qw+I)] 

where 


(29b) 


(30a) 
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=  z*  -  H  (Fi-ii  1  +  Gh»h) 


(30b) 


3.7  Robust  Statistics  &  Stochastic  Approximation 

In  addition  to  simple  tests  described  above,  more  robust  criteria 
can  be  used  to  estimate  variance  and  correlation  coefficients  so  as 
to  make  the  results  less  insensitive  to  outliers  and  have  smaller 
bias  and  variations  [21].  Examples  include  the  bi-weight  method 
[22]  and  the  square  of  the  absolute  median  deviation  for  robust 
variance  estimation  [9,  22]  and  the  rank  correlation  method  for 
estimating  correlation  coefficients  [9]. 

Stochastic  approximation  methods  (SA)  provide  an  online 
framework  for  robust  non-parametric  parameter  estimation  [13], 
system  identification  [1],  and  quantile  estimation  among  others 
[25].  In  particular,  the  Robbins  Monroe  SA  and  the  robustized 
minimum  variance  least-square  method  can  be  used  to  estimate 
the  quantiles  of  the  innovations  sequence.  The  quantiles  can  be 
directly  related  to  the  known  properties  of  the  Gaussian 
distribution.  In  conjunction  with  skewness  and  kurtosis,  it  allows 
characterizing  the  deviation  of  the  data  sequence  from  Gaussian. 

Computation  of  quantiles.  It  can  be  shown  that  to  estimate 
quantiles  of  a  distribution,  the  SA  recursion  is  of  the  form  [14, 
17]: 


(31a) 


where 


S,„  (xp ,  K )  =  —  2sgn(x  -  Yi+m  „ )  +  (1  ■ -  2At )  (3  lb) 

m  "T 

where  Y  is  the  data  sequence,  assumed  to  be  i.i.d,  Ap  is  an 
adaptive  gain  matrix  (based  on  a  nonparametric  estimator)  and 
Ak  is  the  selected  quantile  to  be  estimated.  Note  that  k  quantiles 
can  be  simultaneously  estimated  by  this  procedure.  Further  note 
that  the  computation  of  the  adaptive  Ap  requires  batch  processing 
and  ranking  incurring  a  small  delay.  Alternatively,  a  non-adaptive 
gain  coefficient  could  be  used  without  ranking.  However,  while 
streamlining  processing  it  would  not  possess  the  robust  properties 
of  the  adaptive  SA. 

On-line  density  estimation.  For  recursive  density  estimation,  the 
SA  procedure  is  modified  to  estimate  parameters  of  a  non- 
orthogonal  basis  function  (translated  Beta  density  function) 
approximation  of  an  arbitrary  density  function  [14,  17],  The 
recursion  in  this  case  becomes: 

Vi  =xp--^Ap(xp-D-'</>(Wp))  (32) 

where  t/fw)  is  the  translated  Beta  density  function,  and  w  are  the 
i.i.d  random  samples  from  the  distribution,  and  D  is  band  matrix 
with  2k  +1  nonzero  diagonals  given  by  dik=  <p>. 

This  method  minimizes  the  integral  square  error  between  the 
“model”  and  the  actual  density  to  be  estimated.  A  salient  feature 
of  this  method  is  that  it  does  not  require  storage  or  ranking  and 
can  be  used  directly  on-line  to  estimate  the  “empirical”  density 
function  as  function  of  time.  The  characteristics  of  the  empirical 
density  can  be  compared  with  respect  to  the  normal  density  with 
deviations  computed  by  methods  as  described  above. 

3.8  Illustrative  Examples 

The  same  simulation  scenarios  as  in  Section  2  are  used  here  to 
illustrate  optimality  monitoring  where  the  sliding  a  window 
length  of  40  sec  is  used,  which  is  roughly  equal  to  the  filter 


transient  interval.  Fig.  16  shows  the  test  statistic  of  the  zero  mean 
condition  for  the  innovation  with  upper  and  lower  significant 
bands  on  the  same  plot.  For  a  sample  run,  the  condition  is 
satisfied  over  the  quiescent  periods  without  maneuver.  However, 
it  is  violated  in  all  maneuver  periods. 

Fig.  17  shows  the  skewness  (blue)  and  kurtosis  (green)  estimated 
over  each  sliding  window  of  innovations.  After  a  transient,  both 
estimates  go  down  to  small  values.  Right  after  a  maneuver 
(initiation  and  termination),  large  values  appear,  indicating 
deviation  from  Gaussian.  The  signs  of  skewness  are  consistent  of 
the  directions  of  acceleration-induced  biases.  The  kurtosis  shows 
large  concentration  (value  peaks)  right  after  the  initiation  and 
termination  of  a  maneuver. 

For  a  particular  run.  Figs.  18  and  19  are  the  histogram  of  the 
innovations  in  the  quiescent  and  maneuver  modes,  respectively. 

Fig.  20  shows  the  y2  test  of  the  innovation,  which  satisfies  the 
zero-mean  white  Gaussian  condition  in  quiescent  modes  but 
violates  it  in  maneuver  modes. 

Fig.  21  shows  the  covariance  function.  For  this  particular  data 
window,  except  for  two  small  spikes,  the  covariance  values  at 
other  lags  are  within  the  thresholds,  indicating  no  correlation, 
thus  white  practically. 

The  same  simulation  was  run  for  the  sensor  noise  only  case 
where  the  random  acceleration  noise  is  set  to  zero  v,  =  0  with 
only  deterministic  acceleration  left.  Fig.  22  shows  the  test 
statistic  of  the  zero  mean  condition  for  innovation  with  upper  and 
lower  significant  bands.  Fig.  23  shows  the  skewness  and  kurtosis 
estimates.  Fig.  24  shows  the  autocorrelation  function  with  upper 
and  lower  significant  bands. 

The  results  are  “cleaner”  than  the  previous  case  with  process 
noise  and  all  the  conclusions  about  the  tests  hold. 

Additional  results  for  probability  of  detection  and  false  alarm  rate 
as  well  as  latency  of  these  and  other  test  statistics  will  be 
examined  in  subsequent  papers. 

4  Adaptive  Fusion  &  Management 

In  run  time,  a  tracking  filter  needs  to  be  prepared  to  handle 
possible  changes  in  sensor  error  characteristics  and  target 
maneuvers,  both  leading  to  model  mismatches  and  if  not  dealt 
with  properly,  would  degrade  the  overall  tracking  performance. 

The  optimality  self  online  monitoring  (OSOM)  procedure 
described  in  Section  3  provides  means  to  verify  design 
assumptions  and  validate  operating  conditions.  Any  deviation 
from  the  optimality  prompts  actions.  This  may  range  from 
individual  filter  tuning,  adaptive  sensor  fusion,  active  sensor 
scheduling  and  management,  to  network  coordination. 

4.1  Individual  Filter  Tuning 

It  is  obvious  that  when  the  actual  measurement  error  covariance 
R  is  found  to  differ  significantly  from  the  assumed  value  in  the 
filter,  the  filter  ought  to  be  adjusted  to  reflect  the  new  reality  for 
better  performance. 

The  design  procedure  described  in  Section  2  factors  in  maximum 
possible  maneuver  but  in  a  conservative  manner.  In  contrast, 
proactive  methods  include  reactive  adaptation  and  multiple 
model  estimation  [3],  The  former  includes  single  filter 
adaptation,  variable  dimension  filter,  two-stage  bias  estimation, 
and  input  estimation.  The  latter  is  exemplified  by  the  Multiple 
Model  Adaptive  Estimator  (MMAE)  [20]  and  Interacting 
Multiple  Model  (IMM)  algorithm  [5]. 
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4.2  Adaptive  Fusion  with  Optimality  Monitoring 

These  maneuver  adaptive  filtering  methods  still  experience  biases 
or  lags  in  position  and  velocity  estimates  after  the  initiation  and 
termination  of  a  maneuver,  although  reduced  considerably  as 
compared  to  non-maneuver  filters.  Such  errors  are  primarily 
caused  by  transient  behaviors.  If  these  estimates  are  used  in  track 
fusion,  the  fused  track  may  be  “derailed”  if  different  biases  or 
lags  from  different  sensors  are  not  recognized.  Since  the 
optimality  conditions  no  longer  hold  in  the  innovation,  statistical 
tests  similar  to  those  for  individual  filters  as  described  in  Section 
3  may  be  developed  to  conduct  fuser  autonomous  integrity 
monitoring  (FAIM)  to  ensure  consistent  fusion. 

4.3  Sensor  Management  &  Network  Coordination 

One  way  to  reduce  transient  errors  is  to  obtain  a  direct  estimate  of 
the  target  maneuver  and  to  correct  the  tracking  filtering 
accordingly.  In  the  past,  change  in  target  orientation  from  a 
sequence  of  visual  images  was  used  to  deduce  target  maneuver 
[27].  Recently,  range-Doppler  images  of  a  target  from  a  high 
range  resolution  (HRR)  radar  are  used  to  extract  target  maneuver 
information  [28].  This  involves  change  of  sensor  modes  and 
coordination  of  multiple  sensors  from  the  same  or  different 
platforms. 

As  described  in  this  paper,  the  OSOM  procedure  is  a  bottom-up 
approach  from  individual  filters  through  fusion  centers  to 
network  managers,  which  will  play  a  critical  role  for  quality 
insurance  in  network-centric  layered  sensing. 

5  Conclusions 

In  this  paper,  we  reported  our  initial  study  of  online  monitoring 
of  tracking  optimality  as  means  for  performance  evaluation  and 
adaptation.  The  steady-state  Kalman  filter  in  one  coordinate,  i.e., 
the  a-/5  filter,  was  used  as  an  example  for  illustration.  As  shown 
by  simulation,  the  tracking  performance  could  be  well  predicted 
when  the  models  were  perfectly  known.  However,  it  failed  when 
the  filter  models  did  no  match  reality.  Various  optimality  tests 
applied  to  the  innovation  sequence  could  detect  such  conditions. 

The  target  state  estimate  and  the  corresponding  estimation  error 
covariance  should  be  accompanied  by  an  optimality  indicator 
(consistency  or  integrity  indicator)  when  it  is  used  for  track 
fusion,  sensor  management,  targeting,  and  control.  Our  ongoing 
research  in  this  direction  seeks  for  fast  and  reliable  indicators  and 
metrics  and  their  incorporation  to  target  tracking,  sensor  fusion, 
and  resource  management  algorithms. 
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